00_census_2020_exploratory
From NYC Planning:
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.1.1## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.5 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.0.2 v forcats 0.5.1## Warning: package 'ggplot2' was built under R version 4.1.1## Warning: package 'tibble' was built under R version 4.1.1## Warning: package 'tidyr' was built under R version 4.1.1## Warning: package 'readr' was built under R version 4.1.1## Warning: package 'purrr' was built under R version 4.1.1## Warning: package 'dplyr' was built under R version 4.1.1## Warning: package 'stringr' was built under R version 4.1.1## Warning: package 'forcats' was built under R version 4.1.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(tidycensus)## Warning: package 'tidycensus' was built under R version 4.1.1library(jsonlite)## Warning: package 'jsonlite' was built under R version 4.1.1##
## Attaching package: 'jsonlite'## The following object is masked from 'package:purrr':
##
## flattenlibrary(magrittr)## Warning: package 'magrittr' was built under R version 4.1.1##
## Attaching package: 'magrittr'## The following object is masked from 'package:purrr':
##
## set_names## The following object is masked from 'package:tidyr':
##
## extractlibrary(sf)## Warning: package 'sf' was built under R version 4.1.1## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1library(tigris)## Warning: package 'tigris' was built under R version 4.1.1## To enable
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.##
## Attaching package: 'tigris'## The following object is masked from 'package:tidycensus':
##
## fips_codeslibrary(censusapi)## Warning: package 'censusapi' was built under R version 4.1.1##
## Attaching package: 'censusapi'## The following object is masked from 'package:methods':
##
## getFunctionreadRenviron("~/.Renviron")
get_decennial(
geography = "tract",
state = "NY",
variables = "P1_001N",
year = 2020,
geometry = T
)## Getting data from the 2020 decennial Census## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.## Using the PL 94-171 Redistricting Data summary file##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|======================= | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|=================================== | 49%
|
|=================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|=============================================================== | 89%
|
|=============================================================== | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|===================================================================== | 99%
|
|======================================================================| 100%## Note: 2020 decennial Census data use differential privacy, a technique that
## introduces errors into data to preserve respondent confidentiality.
## i Small counts should be interpreted with caution.
## i See https://www.census.gov/library/fact-sheets/2021/protecting-the-confidentiality-of-the-2020-census-redistricting-data.html for additional guidance.
## This message is displayed once per session.#load_variables(2020, "pl")
varlist <- fromJSON("https://api.census.gov/data/2020/dec/pl/variables.json", simplifyDataFrame = T) %>%
extract2("variables") %>%
as.data.frame() %>%
mutate_all(~as.character(.)) %>%
pivot_longer(everything(), names_sep = "\\.",
names_to = c("varname", "type"))Load census data
ny_t_2020 <- st_read(dsn = "C:/Users/jg6849/Documents/Github/census-2020/data/tl_2020_36_tract/tl_2020_36_tract.shp") %>%
filter(COUNTYFP %in% c("005", "047", "061", "081", "085"))## Reading layer `tl_2020_36_tract' from data source
## `C:\Users\jg6849\Documents\Github\census-2020\data\tl_2020_36_tract\tl_2020_36_tract.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 5411 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -79.76259 ymin: 40.47658 xmax: -71.77749 ymax: 45.01586
## Geodetic CRS: NAD83ny_t_2020 ny_t_2010 <- st_read(dsn = "C:/Users/jg6849/Documents/Github/census-2020/data/tl_2010_36_tract10/tl_2010_36_tract10.shp") %>%
filter(COUNTYFP10 %in% c("005", "047", "061", "081", "085")) %>%
rename(GEOID = GEOID10)## Reading layer `tl_2010_36_tract10' from data source
## `C:\Users\jg6849\Documents\Github\census-2020\data\tl_2010_36_tract10\tl_2010_36_tract10.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 4919 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -79.76259 ymin: 40.4774 xmax: -71.77749 ymax: 45.01586
## Geodetic CRS: NAD83ny_t_2010#ny_t_2020 <- tracts("NY", county = c("New York", "Kings", "Queens", "Richmond", "Bronx"),
# year = 2020)
#
# ny_t_2010 <- tracts("NY", county = c("New York", "Kings", "Queens", "Richmond", "Bronx"),
# year = 2010)
# ny_t_2020
# ny_t_2010ny_t_2020_filt <- ny_t_2020 %>%
filter(AWATER < ALAND)
ny_t_2010_filt <- ny_t_2010 %>%
filter(AWATER10 < ALAND10)
tract_intersection <- function(gdf1, gdf2, id) {
gdf1 %>%
filter(GEOID == id) %>%
st_intersection(filter(gdf2, GEOID == id))
}
tract_difference <- function(gdf1, gdf2, id) {
gdf1 %>%
filter(GEOID == id) %>%
st_difference(filter(gdf2, GEOID == id))
}
ny_t_int <- map_dfr(unique(ny_t_2020_filt$"GEOID"), ~tract_intersection(ny_t_2020_filt, ny_t_2010_filt, .))
ny_t_intny_t_diff_2020 <- map_dfr(unique(ny_t_2020_filt$"GEOID"), ~tract_difference(ny_t_2020_filt, ny_t_2010_filt, .))
ny_t_diff_2020ny_t_diff_2010 <- map_dfr(unique(ny_t_2010_filt$"GEOID"), ~tract_difference(ny_t_2010_filt, ny_t_2020_filt, .))
ny_t_diff_2010stopifnot(st_crs(ny_t_2020_filt) == st_crs(ny_t_2020_filt))Intersection of census tracts between 2020 and 2010
ggplot(data = ny_t_int) +
geom_sf(fill = "gray", lwd = 0.5) +
theme_classic() +
theme(
axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank()
)
Difference in census tracts between 2020 and 2010
ggplot() +
geom_sf(data = ny_t_diff_2020, fill = "blue", color = "blue", alpha = 0.5) +
geom_sf(data = ny_t_diff_2010, fill = "green", color = "green", alpha = 0.5) +
theme_classic() +
theme(
axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank()
)
ny_t_diff_comp <- ny_t_diff_2020 %>%
mutate(area_2020diff = st_area(geometry)) %>%
as.data.frame() %>%
select(c("GEOID", "area_2020diff")) %>%
full_join(ny_t_diff_2010 %>%
mutate(area_2010diff = st_area(geometry)) %>%
as.data.frame() %>%
select(c("GEOID", "area_2010diff")), by = "GEOID") %>%
left_join(ny_t_2020_filt, ., by = "GEOID") %>%
mutate(area_2020_total = as.numeric(st_area(geometry))) %>%
group_by(GEOID) %>%
mutate(area_ratio = as.numeric((area_2010diff + area_2020diff) / area_2020_total)) %>%
arrange(desc(area_ratio)) %>%
select("area_ratio", everything()) %>%
ungroup() %>%
mutate(area_ratio = case_when(is.na(area_ratio) ~ 0,
TRUE ~ area_ratio))
ny_t_diff_compstopifnot(eeptools::isid(ny_t_diff_comp, "GEOID"))Scored difference
ggplot() +
geom_sf(data = ny_t_diff_comp, aes(fill = area_ratio)) +
theme_classic() +
theme(
axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank()
) +
scale_fill_continuous(high = "#132B43", low = "#56B1F7")
library(leaflet)## Warning: package 'leaflet' was built under R version 4.1.1library(mapview)## Warning: package 'mapview' was built under R version 4.1.1#m <- mapview(ny_t_diff_comp)
#m
ny_t_diff_comp <- ny_t_diff_comp %>%
filter(!is.na(area_ratio)) %>%
st_transform(crs = st_crs(4326)) %>%
filter(area_ratio > 0)
pal_rev <- colorNumeric(colorRamp(c("#40899A", "#FFFFFF"),
interpolate="spline"),
domain = ny_t_diff_comp$area_ratio)
pal <- colorNumeric(colorRamp(c("#40899A", "#FFFFFF"),
interpolate="spline"),
domain = ny_t_diff_comp$area_ratio, reverse = T)
pal## function (x)
## {
## if (length(x) == 0 || all(is.na(x))) {
## return(pf(x))
## }
## if (is.null(rng))
## rng <- range(x, na.rm = TRUE)
## rescaled <- scales::rescale(x, from = rng)
## if (any(rescaled < 0 | rescaled > 1, na.rm = TRUE))
## warning("Some values were outside the color scale and will be treated as NA")
## if (reverse) {
## rescaled <- 1 - rescaled
## }
## pf(rescaled)
## }
## <bytecode: 0x000000000c6281f8>
## <environment: 0x0000000004c11240>
## attr(,"colorType")
## [1] "numeric"
## attr(,"colorArgs")
## attr(,"colorArgs")$na.color
## [1] "#808080"leaflet(data = ny_t_diff_comp) %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolygons(data = ny_t_2010_filt,
color = "black", stroke = T, weight = 0.2,
fill = T, fillColor = "grey",
group = "Census Tracts 2010") %>%
addPolygons(color = "gray", stroke = T,weight = 0.2,
fill = T, fillColor = ~pal(area_ratio), fillOpacity = 1,
group = "Census Tracts 2020") %>%
addLegend(position = "topright",
pal = pal_rev,
values = ~area_ratio,
title = "Area ratio",
labFormat = labelFormat(transform = function(area_ratio)
sort(area_ratio, decreasing = TRUE))) %>%
addLayersControl(overlayGroups = c("Census Tracts 2010", "Census Tracts 2020"))## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'